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Central to the functioning of a living cell is its ability to control the readout or expression of 
information encoded in the genome. In many cases, a single transcription factor protein activates or 
represses the expression of many genes. As the concentration of the transcription factor varies, the 
target genes thus undergo correlated changes, and this redundancy limits the ability of the cell to 
transmit information about input signals. We explore how interactions among the target genes can 
reduce this redundancy and optimize information transmission. Our discussion builds on recent work 
[Tkacik et al, Phys Rev E 80, 031920 (2009)], and there are connections to much earlier work on the 
role of lateral inhibition in enhancing the efficiency of information transmission in neural circuits; 
for simplicity we consider here the case where the interactions have a feed forward structure, with 
no loops. Even with this limitation, the networks that optimize information transmission have a 
structure reminiscent of the networks found in real biological systems. 



I. INTRODUCTION 

The genomes of even the smallest bacteria encode the 
structure of hundreds of proteins; for complicated organ- 
isms like us, the number of different proteins reaches into 
the tens of thousands [T]. During the course of its life, 
each cell has to control how many copies of each protein 
molecule are synthesized [2J. These decision about the 
expression of genetic information occur on many scales, 
from (more or less) irreversible decisions during differ- 
entiation into different tissue types down to continuous 
modulations of the number of enzyme molecules in re- 
sponse to varying metabolic needs and resources [3]. As 
with many biological processes, the control of gene ex- 
pression depends critically on the transmission of infor- 
mation. Because the relevant information is represented 
inside the cell by signaling molecules at low concentra- 
tions, the irreducibly stochastic behavior of individual 
molecules sets a physical limit to information transmis- 
sion in |5j [6l [7] . A more precise discussion is motivated 
by the emergence of experiments which measure the noise 
in gene expression [HJ|ll[Tni[IU[Tl[ni[Tl]. In this paper, 
building on our previous work |15j . we explore how cells 
can structure their genetic control circuitry to optimize 
information transmission in the presence of these physi- 
cal constraints, thus maximizing the control power that 
they can achieve while using only a limited number of 
molecules. 

Although the problems of information transmission in 
genetic control are general, it is useful to have a con- 
crete example in mind. In developing embryos, infor- 
mation about position — and hence, ultimately, fate in 
the developed organism — is encoded by spatially varying 
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concentrations of 'morphogen' molecules [inillT]. These 
molecules often are transcription factors, and positional 
information then is transmitted to the expression levels 
of the target genes [TTl [THl HH] • In this scheme, maximiz- 
ing information transmission maximizes the richness of 
the spatial patterns which the embryo can construct. 

In recent work we have considered the problem of max- 
imizing information transmission from a single transcrip- 
tion factor to multiple target genes, in the case where the 
targets are non-interacting [15]. Already this problem 
generates some structures that remind us of real genetic 
control networks. But, when a single transcription factor 
controls many genes independently, the signals carried 
by those genes necessarily are redundant [20]. To fully 
optimize the information which can be carried by a given 
number of molecules, there must be interactions among 
the target genes which reduce this redundancy. Our goal 
in this paper is to derive the form of these interactions. 
We emphasize that this is an ambitious project: we are 
trying to find the topology and parameters of a genetic 
network from first principles, constrained only by the lim- 
ited concentration of all the relevant molecules. To sim- 
plify our task, we start here with the case in which the 
network of interactions has no closed loops, and return 
to the full problem in subsequent papers. 

When we search possible networks for the ones that op- 
timize information transmission, it is convenient to break 
this search into two parts. First, we enumerate network 
topologies, and then we search parameters within each 
topology. It is conventional to encode in the network 
topology the sign of the interactions, so that changing 
the 'shape of an arrow' from activation to repression is 
counted as a qualitatively different network (see Fig. [TJ . 
Even if we assume that each protein can act either as 
activator or repressor for all its targets [H], with one 
transcription factor controlling K interacting genes in 
a feed forward network, there are 2^ possible network 
topologies. For each of these possibilities we find a single 
well defined optimum for all the continuously adjustable 



2 



ACTIVATION: REPRESSION: 




A > B A 1 B 

b. C 




Qi ' 92 

FIG. 1: Schematic models of transcriptional regulation, (a.) 
Transcription factors proteins (A) can either activate (right) 
the expression of genes (B) by recruiting the RNA poly- 
merase, or repress (left) the expression of genes by block- 
ing the polymerase. Throughout the paper we will use the 
diagrammatic representation, in which an arrow depicts acti- 
vation and a blunted arrow depicts repression, (b.) A feed 
forward network in which the input c regulates the expression 
levels of gi and (72, and gi also regulates 32; the structure is 
feed forward because g2 does not feedback on gi. 



parameters in the system, and these optimal parameter 
values are determined only by the number of available 
molecules. Thus, the structures of these networks really 
are derivable, quantitatively, from first principles. For re- 
lated approaches to information flow in biochemical and 
genetic networks, see Refs [231 [Ml 123 US 127] . 

Among the possible network topologies, there are some 
in which interaction strengths are driven to zero by the 
optimization. Among the nontrivial solutions, however, 
we find that the different topologies achieve very similar 
information capacities. Thus, it seems that there is not a 
single optimal network, but rather an exponentially large 
number of nearly degenerate local optima. All of these 
local optima are networks with competing interactions, 
so that a single target gene is activated and repressed 
by different inputs. Phenomenologically, if we view the 
expression level of each target gene as a function of the 
input transcription factor concentration, these compet- 
ing interactions lead to non-monotonic input/output re- 
lations. It is this non-monotonicity which allows the col- 
lection of target genes to explore more fully the space of 
expression levels, thus enhancing the capacity to transmit 
information. Although we consider only networks with- 
out feedback, these non-monotonic responses already are 
reminiscent of the patterns seen in real genetic networks. 



II. SETTING UP THE PROBLEM 

The optimization problem that we are addressing in 
this paper is a generalization of the problem discussed 
in our previous work [15j . To make this paper self- 
contained, however, we begin by reviewing the general 
formulation of the problem. We are interested in sit- 
uations where a single transcription factor controls the 
expression level of several genes, labelled i = 1,2, - ■ ■ ,K. 
The expression levels of these genes, gi, carry informa- 
tion about the concentration c of the input transcription 
factor. It is this information /(c; {^i}) that we suggest is 
optimized by real networks. 

To derive the consequences of our optimization princi- 
ple, we require several ingredients: we need to describe 
the space of possible networks, we need to relate /(c; {^i}) 
to some calculable properties of these networks, we need 
to put these pieces together to form the objective func- 
tion for our optimization problem, and finally we need to 
perform the optimization itself. In the interests of pro- 
ceeding analytically as far as possible down this path, 
before resorting to numerics, we adopt several approx- 
imations: we consider the case where expression levels 
have reached steady state values consistent with the in- 
put transcription factor concentration, we assume that 
the noise in the system is small, and we restrict our at- 
tention to networks that have a feed forward structure. 
We have discussed the first two approximations in our 
previous work [15j, and our focus here on feed forward 
networks is intended as a useful intermediate step be- 
tween the case of non-interacting target genes (as in Ref 
[T5|) and the case of arbitrary interactions, to which we 
will return in a subsequent paper. 

A. Describing the regulatory interactions 

The expression levels of genes are determined by a 
complex sequence of events, especially in eukaryotic cells 
where even the transcription of DNA into messenger 
RNA involves a complex of more than one hundred pro- 
teins [22]. As in much previous work, we abstract from 
this complexity to say that expression levels are set by 
a balance of synthesis and degradation reactions; in the 
simplest case only synthesis is regulated, and degradation 
is a first order process. Then the deterministic kinetic 
equations that the govern the expression levels are 

rlC- 1 

w/i(c;{Gj})~-Gi, (1) 

where r is the lifetime of the gene products against degra- 
dation, Tmax is the maximum number of molecules per 
second that can be synthesized, and /i(c;{Gj}) is the 
"regulation function" that describes how the synthesis 
rate is modulated by the other molecules in system. The 
regulation functions, which are positive and range be- 
tween zero and one, will be different for every gene, but 
we assume that the lifetimes and maximum synthesis 
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rates are the same for all of the genes [55] . In this formu- 
lation, Gi counts the number of molecules of the protein 
coded by gene i, while it is conventional to define c as 
a concentration. Before we are done it will be useful to 
normalize these quantities. 

The steady state expression levels, Gi{c), are the solu- 
tions of the simultaneous nonlinear equations 



the spectral density of synthesis noise, so we have 



Gi(c) = wr/i(c;{Gj(c)}). 



(2) 



In general there can be multiple solutions, corresponding 
to a network that has more than one stationary state; 
as noted above, we will confine our attention here to 
networks with a feed forward architecture, in which this 
can't happen. In Figjlja we present an example of a feed 
forward network, where an input c regulates the expres- 
sion of genes gi and 172 • The product proteins of 171 also 
regulate gene (72, but (72 does not regulate 171. Formally, 
the 'no loops' condition defining feed forward networks 
means that there is some assignment of the labels i so 
that /i only depends upon Gj for j < i, and hence, as will 
be important below, the matrix dfi/dGj is lower trian- 
gular. 

A deterministic system with continuous inputs and 
outputs can transmit an infinite amount of information. 
In reality, information transmission is limited by the 
discreteness and randomness of the individual molecu- 
lar events. If this noise is small, it can be described 
by adding Langevin 'forces' to the deterministic kinetic 
equations [H HUSH EH ES] • To the extent that synthe- 
sis and degradation rates can be written as depending 
on the instantaneous concentrations or expression lev- 
els, as we have done in Eq ([T|), then the system has no 
'hidden' memory, and the Langevin terms which describe 
the noise in the system should have effectively zero cor- 
relation time. We also will assume that the noise in the 
synthesis and degradation of different gene products are 
independent. Then we can write 



dGi 



1 



(3) 
(4) 



where the A/i are the spectral densities of the noise forces. 

The spectral densities of noise have several components 
which add together. First, if we take Eq (jlj seriously, it 
describes synthesis and degradation as simple first order 
processes; the terms which appear in the deterministic 
equation as rates should then be interpreted as the mean 
rates of Poisson processes that describe the individual 
molecular transitions. Then we have 



AA/') = w/i(c;{Gj}) + -Gi. 

T 



1 



(5) 



We are interested in small fluctuations around the steady 
state, so in this expression we should set Gj = Gi(c). 
Further, we can account for the possibility that synthesis 
is a multi-step process by adding a 'Fano factor' F to 



A/;' 



(1) 



i^w/i(c;{Gj(c)})-H-Gi(c) (6) 

T 

l + F 



Gi(c). 



(7) 



The Fano factor is a measure of the dispersion of the 
probability distribution, the ratio of the variance to the 
mean; here the relevant random variable is the number 
of reaction events in a small window of time. As an 
example, if synthesis occurs in bursts, then we expect 
F > 1 [Sni [311 135]. Burst of both proteins and niRNAs 
have now been experimentally observed in a number of 
bacterial and eukaryotic systems [51 [351 [SZl [3S] . 

A second source of noise comes from the diffusive ar- 
rival of transcription factor molecules at their targets. To 
describe this, we should remember that the regulation of 
transcription depends on events that occur in a very small 
volume, of linear dimension £. The relevant concentra- 
tions are those in this volume, or more colloquially at the 
point where the transcription factors bind; binding and 
unbinding events acts as localized sources and sinks for 
the diffusion equation which describes the spatiotempo- 
ral variations in concentration. Intuitively, this coupling 
generates noise because the diffusion of the transcription 
factors into the relevant volume is a stochastic process, 
so that even with the global concentration fixed there 
are local concentration fluctuations [35^ . Analysis of this 
problem [4][5l|40] shows that, so long as the time for dif- 
fusion through a distance £ is short, the net effect of the 
diffusive fluctuations is captured by a Langevin term 



(2) 



dfijc; {Gj}) 
dc 



c 



(8) 



where D is the diffusion constant; again wc should eval- 
uate this at Gj = Gj (c) to describe the small fluctuations 
around the steady state. 

The same arguments which apply to the input tran- 
scription factor also apply to each of the target gene prod- 
ucts when they act as transcription factors. If there are 
Gi molecules of gene product i, then these proteins are 
present at concentration Gi/il, where is the relevant 
volume [3T] , and the analog of Eq dsl is 
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a/i(c;{Gk})V (Gk/f^) 



diG^/n) 

dMc; {Gk}) 
dG. 
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2 

Gk. 



(9) 
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Finally, in putting all these terms together, it is con- 
venient to measure expression levels in units of the 
maximum mean expression level, which is Gj — 7',„ax''" 
molecules; that is, we define gi = Gi/(rmaxT). Then we 
have 
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T— — 

dt 



{.9j})-5i+6(i), 



(5ij(5(t - t')T 



C 



E 



{ffk}) 
9gk 



.9k 



(11) 

(12) 



{ffk=gk(c)} 
(13) 



The parameter combinations which appear in these 
equations have simple interpretations. If _F = 1, the syn- 
thesis and degradation of each gene product molecule re- 
ally is an independent Poisson process, and hence r^^y^r 
is not just the maximum number of molecules that can 
be made (on average), but also the maximum number 
of independent molecules, Ng in the notation of Ref 
[15j . For more complex processes, where F ^ 1, we can 
write Ng = 2ri„ax''"/(l + F). Similarly, the combination 



r-aia-xT /U, is a concentration, the maximum (mean) con- 
centration of the gene products. Since these are acting 
as transcription factors, we will assume that this maxi- 
mum concentration is also the maximum concentration 
of the input transcription factor, Cmax- As discussed pre- 
viously [15 , there is a natural scale of concentration, 
Co = Ng/^Dir), and if we use units in which cq = 1 
we can write 



3i(c) 



/ 9/i(c;{gk}) 
I dc 



( dfijc; {gk}) 



9k 



(14) 



{9k=Sk(c)} 



r 



To complete our description of the system we need to 
specify the regulation functions /i(c; {.9j}). 

B. The regulation functions 

The central event in transcriptional regulation is the 
binding of the transcription factor (s) to their specific sites 
along the DNA. These binding sites must be close enough 
to the start of the coding sequence of the target gene 
that binding can influence the initiation of transcription 
by the RNA polymerase and its associated proteins. In 
the simplest geometric picture, which might actually be 
accurate in bacteria [HI |35] HH US] , we can imagine the 
molecular interactions being so strong that the normal- 
ized mean rate of transcription — what we write as the 
regulation function /i(c;{gj}) — is essentially the proba- 
bility that certain binding sites are occupied or unoc- 
cupied. If n transcription factors bind cooperatively to 
nearby sites, and binding activates transcription, then we 
would expect 

c" 

and if the binding represses transcription, so that the 
mean rate is proportional to the fraction of empty sites. 



we should have 

/(^) - (16) 

In each case, the constant K measures the concentra- 
tion for half-occupancy of the binding sites, and so 
F ~ —ksTlnK can be interpreted as a binding en- 
ergy for the transcription factor to its specific site along 
the genome. In the biological literature these models 
are called Hill functions, after Hill's classical discussion 
of oxygen binding to hemoglobin [46^. While it is con- 
venient to distinguish between activators and repressors, 
we note that, within the Hill function model, we can pass 
smoothly from one to the other by allowing n to change 
sign. 

The Hill function model identifies the transcription of 
a gene as a logical function of binding site occupancy, so 
that the mean rate of transcription will be proportional 
to the probability of the sites being occupied (or, in the 
case of repressors, not occupied). Thus there is a natural 
generalization when multiple transcription factors con- 
verge to control a single gene, as in Fig [2] the rate of 
transcription should be proportional to logical functions 
on the occupancy of multiple binding sites. Thus, if the 
input transcription factor and the gene j converge to ac- 
tivate gene i, transcription could depend on the AND of 
the two binding events; if these are independent, then 



5 



the regulation function becomes 



(17) 



Alternatively, transcription could depend on the OR of 
the two binding site occupancies, in which case 



/i(c,5j) 



(18) 



and we can construct functions that interpolate or mix 
AND and OR. In general we find that the pure AND 
functions transmit the most information. Thus we will 
write 



/i(c,{5j}) 



n 



3i 



(19) 



where the restriction j < i confines our attention to feed 
forward networks with no loops. If ny > 0, 



then Kij 







means that gene j does not regulate gene i; if < then 
the condition for no interaction is Kij — > oo. 

An alternative model for regulation envisions a large 
transcriptional complex that can be in one of two states, 
and the rate of transcription is proportional to the popu- 
lation of the 'active' state; transcription factors act by 
binding and shifting the equilibrium between the two 
states. This is essentially the model proposed by Monod, 
Wyman and Changeaux for allosteric enzymes and co- 
operative binding [47]. The important idea is that, in 
each state, all effector molecules bind independently, but 
the binding energies are different in the two states. By 
detailed balance, this difference in binding energy means 
that binding will shift the equilibrium between the two 
states. In the case of one transcription factor at concen- 
tration c, with n binding sites, the probability of being 
in the active state is given by: 



/(c) 



(1 + c/K,^y 



L"(l + c/i^off)"+(l+c/i^on)"' 



(20) 



where kgT In L is the free energy different between the 
two states in the absence of any bound molecules. If 
binding favors the active state, and the gene has a very 
small probability to be "on" in the absence of input {L ^ 
1), then taking K^n <C K^s reduces Eq (20 1 to: 



/(c) 



(1 + c/KqX 

1 



1 + exp 



-nln 



l + Ci/2/Ka 



(21) 
(22) 



where we have introduced the constant Ci/2 to simplify 
notation, i = (1 + Ci/2/Koa), and Ci/2 is the input con- 
centration for half maximal activation [48]. Note that 
activation and repression differ only in the sign of n. 

When more than one type of protein regulates expres- 
sion (Fig^), each binding event makes its independent 



contribution to shifting the equilibrium between active 
and inactive states. Thus, in the same limit of very un- 
equal binding constants to the two states, we can write 



/i(c,{ffj}) 
Fi{c,{9i}) 



1 



l+exp[-f^i(c,{5j})] 
1 + c/Ki 



(23) 



rii In 



1 + c« /ifi 



j<i 



1 



+ .9j(c« )/i^ij 



.(24) 



For a single transcription factor, there is not much dif- 
ference between the Hill and MWC models: the classes 
of functions are very similar, and when we go through 
the optimization of information transmission, the opti- 
mal parameter values are such that the input /output 



b: 



AND 




FIG. 2: Regulation of transcription in the presence of two 
transcription factors, (a.) Two types of transcription factors, 
squares and triangles, can bind to the regulatory region of the 
gene. This situation does not uniquely determine the form of 
regulation, (b.) The Hill regulatory model with the AND 
logic for the binding of the two types of transcription factors. 
In this case the both types of transcription factors must be 
bound to initiate transcription. The response of the gene is 
proportional to the product of the activation by each tran- 
scription factor separately, (c.) The MWC model for model 
for regulation. The transcriptional complex has two states, 
active and inactive, and binding of any type of transcription 
factor shifts the equilibrium between the two states. 
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relations are nearly indistinguishable. But, as we shall mation that the expression levels carry about the input 
see, when multiple transcription factors converge, the Hill transcription factor concentration. This requires only a 
and MWC models lead to significantly difl^erent behav- modest generalization of the discussion in Ref [TS], but, 
iors. This raises the possibility that the molecular details in the interest of clarity, we review the derivation, 
of the regulatory process are still 'felt' at the more macro- 
scopic, phenomcnological level, but this is a huge set of 
questions that we leave aside for the moment, focusing 
just on these two models. 

In the limit that noise is small, the solution of the 
Langevin equations in Section [ll A| has the form of Gaus- 
sian fluctuations around the steady state. Formally, the 
probability distribution of expression levels given the con- 
centration of the input transcription factor is 



C. Information transmission 

We are interested in connecting the properties of the 
networks described in the previous section to the infor- 



Pi{9i}\c) 



exp 



1 1 ^ 

- logdet(/C) - 7^^{9i- ffi(c)) ^ij (ft - 3j(c)) 



The mean expression level of each gene, 5i(c), defines an 
input/output relation for that gene, and the matrix /C is 
the inverse covariance matrix of the fiuctuations or noise 
in the expression levels at fixed input. 



((.9i - .9i(c)) (.9j 



(26) 



note that /C is a function of c. The fact that the ex- 
pression levels "carry information" about the input tran- 
scription factor concentration means that it should be 
possible to estimate c from knowledge of the {gi}. From 
Bayes' rule, the distribution of concentrations consistent 
with some set of expression levels is given by 



Pi{9i}\c)P{c) 
Pi{9i}) ■ 



(27) 



In the limit that noise is small {K, is large), this too is an 
approximately Gaussian distribution. 



picm) 



1 



exp 



{c-c^{{gi})y 
2a? 



(28) 



i,j = l 



(25) 



r 



where c*({gi}) is the most likely input given the output, 
and the effective variance is determined by 



K 



^gi(c) ^ dgijc) 
dc dc 



(29) 



c=c*({ff,}) 



In general it is difficult to flnd an explicit expression for 
c*({5i}), but we will see that this is not crucial. We can 
use these ingredients to calculate the amount of informa- 
tion, in bits, that the expression levels {^i} provide about 
the input concentration c. 

The mutual information between c and {gi} is defined, 
following Shannon [321 [SD] , as 



/(c;M) = /.c/d-,P(c,M)log, [^^1^" 

(30) 

We can rewrite this as a difference in entropies. 



I{c-{g,})^- j dcP{c)log,P{c)- J d'^gPiigi}) - J dc P{c\{gi}) log, P{c\{gi}) 



(31) 



where the first term is the entropy of the overall distribu- 
tion of input concentrations, and the second term is the 
(mean) entropy of the distribution of inputs conditional 
on the output. With the Gaussian approximation from 



Eq ( 28 1 , we can evaluate 



dcP{c\{g,})log,P{c\{g,}) = -log, [2TTeal{{g,})] . 

(32) 
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Further, when the noise is small, we can replace an av- 
erage over expression levels by an average over inputs, 
setting the expression levels to their mean values along 
this path: 

J d^gPi{g,m--]^ J dcP{c)Y[S{g,-Uc))[---]. 

(33) 

Putting these pieces together, we find an expression for 
the information in the low noise limit, 



He; {5i}) 



dcP{c) log2P(c) 

1 



dcP(c) log2 [27:eaU{Uc)})] ■ 

(34) 



As Equation (34 1 makes clear, the information trans- 
mitted through the regulatory network depends on the 
input/output relations, on the noise level, and on the 
distribution of inputs P(c). The cell can optimize infor- 
mation transmission by adjusting this input distribution 
to match the characteristics of the network [SlIT], but this 
matching is constrained by the cost of making the input 
molecules. We can implement this constraint by trading 
bits against a cost function that counts the mean num- 
ber of molecules, or more simply by insisting that c is 
bounded by some maximum concentration Cmax- In Ref 
[TS] we have shown that these different ways of imple- 
menting the constraint given essentially identical results, 
so we use the Cmax approach here. We still need a con- 
straint to force the normalization of P(c), so we should 
maximize the functional 



£ = /(c, {gi})-A J dcP{c). 



(35) 



The solution to this problem can be written, using the 
expressions for /(c, {^i}) from above, as 



P*(c) cx 



1 



'^c({ffi(c)}) 

K 



(36) 



1 
Z 



1 



27re ^-^ dc 

ij = l 



E 



'^'j^^^ dc 



1/2 



,(37) 



Finally, the information itself, evaluated with the optimal 
input distribution, is 



r(c; {5i})=log2 Z. 



(39) 



This maximal mutual information still depends on all the 
parameters that describe the network, and our goal is to 
find the parameter values that maximize /* . To do this, 
it is useful to be a little more explicit about how the 
parameters enter into the computation of Z. 

D. Putting the pieces together 



To evaluate Z in Eq ( 38 1 , we need to know the inverse 
covariance matrix /Cy that describes the noise in gene 
expression. This can be calculated from the Langevin 
equations in Section |II A[ To do this, we recall that we 
are looking at the small noise limit, so we linearize the 
Langevin equations around the steady state gi — gi{c), 
and then Fourier transform. With gi{t) ~ gi + dgi{t), and 



Sgiit) 



duj 
— f 
27r 



*%(c^), 



(40) 



Eq (111 becomes 



l-iujT 




' Sgi{^) ' 




"li(^)" 


— 021 1 — iujT 




6g2{uj) 






-031 -032 1 - iujT .. 




(5.93 (w) 







where 



, dMc,{gi,}) 

'ij = a 



(41) 



(42) 



{Sk=Sk(c)} 



In this expression we explicitly restrict ourselves to net- 
works with a feedforward architecture, which (as noted 
above) means that we can assign labels i so that gene j 
regulates gene i only if j < i. We recall from Eq ( 14) that 



where the normalization constant 



dc 



K 

E 



^5i(c) ^ /„^c^5j(c) 



27re ^—^ dc 



-/Cij(c)- 



dc 



1/2 



(6(i)ej(0) = '5(i-t')A^ij, 



(38) where the noise matrix is diagonal. 



(43) 



f dMc;{g^}) 
\ dc 



k 



( dMc;{g^}) 



5k 



(44) 



{9k=Sk(c)} 
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The delta function in time means that the noise is white, 
{Uuj)i;{u;')) = 2TTd{oj - c^OMj. (45) 



Equations (41 1 are of the form 

i(c.)-<5i(a;) = e>), 



(46) 



where A^uj) is the matrix describing the linearized dy- 
namics, Sg{uj) — {(5gi((jj), 5g2(cj), • • •}, and similarly for 
f((jj). Evidently Sg{uj) = A^^{uj)£^{u!). We are interested 
in calculating the covariance matrix of the fluctuations 



Sg. Since these are stationary, in the frequency domain 
we will have 

{SU^)S~g;{u:')) = 2ttS{oj - c^')^ij(^), (47) 

and the equal time correlations that we want to calculate 
are related to the power spectrum through 



duj 
2ti 



S'ij(w). 



(48) 



Following through the algebra, we have 



[A-\u^%M-\^%n.{U^)€.{^')) 
[A-\u)],^[A-\oj')]*^^N^^2^5{u: - u:') 



[A-i(a;)]ikiVk,n[A-i(c.)];„,= 
I 



A 



(49) 
(50) 

(51) 



where denotes the Hermitian conjugate of the matrix 



M 



The covariance matrix 
Eq ( 26 1 , then takes the form 



which is the inverse of /Cy in 



dij 
2^ 



A 



(52) 



To complete the calculation of Z we need to do this fre- 
quency integral, and then an integral over the concentra- 
tion c, as in Eq ( 38 1 . Some details of these calculations 
are given in the Appendix. Here we draw attention to 
two key points. 

First, because the noise ^ is white — that is, TV is in- 
dependent of frequency — all of the structure in the fre- 
quency integral comes from the inverse of the matrix 
A{lo). Since the matrix elements of A are linear in fre- 
quency, this results in an integrand which is a rational 
function of w, so the integrals can all be done by closing 
a contour in the complex lo plane. Again, examples are 
in the Appendix. 

Second, and most importantly, when we finish com- 
puting Z it depends on all the parameters we want to 
optimize — e.g., in the Hill description of the regulation 
functions, all of the half-maximal points K-^ and the co- 
operativities n-^ — and in addition it depends on the maxi- 
mal concentration of the input transcription factor, Cmax- 
But there are no other parameters in the problem. Thus, 
if we search for the maximum of Z , and hence the opti- 
mum information transmission, we will find the parame- 
ters of the network as a function of the available number 
of molecules, with no remaining arbitrariness. 



III. OPTIMAL NETWORKS AND THEIR 
PARAMETERS 

In previous work [IB], we considered the case of many 
noninteracting genes regulated by a single transcription 
factor. Optimizing information transmission in this sys- 
tems leads to two qualitatively different families of solu- 
tions, depending on the available dynamic range Cmax/co- 
At small Cmax, the optimal strategy is for the different 
output genes to respond in the same way, allowing the 
system to make multiple independent "measurements" of 
the input concentration; in this case the target genes are 
completely redundant. At large Cmax, a tiling solution 
appears, in which the genes are turned on sequentially at 
different concentrations; the precise setting of the thresh- 
olds is determined by a compromise between minimizing 
the noise and maximizing the use of the dynamic range. 
Even in this regime, however, the target genes are (par- 
tially) redundant. Here we explore how the interactions 
between target genes can be tuned to reduce redundancy 
and increase the capacity of the system to transmit in- 
formation. 



A. Two target genes 

We first consider feed forward networks in which an 
input transcription factor at concentration c regulates 
the expression of two genes, with normalized expression 
levels gi and g2; in addition, the product proteins of gene 
1 regulates the expression of gene 2. Within this class of 
networks, we can calculate the information transmission 
following the outline above, and then search for optimal 
values of all the parameters. 
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FIG. 3: Locally optimal networks for fou r topologies of the two gene feed forward network s, w ith Cmax/co ~ 10. Top panels 
are for the case of the Hill model, Eq (19l, and bottom panels are for the MWC model, Eq (24l. As explained in the text, the 
optimal parameters for the A1A2 — A12 and the R1R2 — A12 networks correspond to zero interaction between the target genes 
so these solutions also provide a reference for input/output relations in the optimal non-interacting networks. Note that the 
non-interacting solutions are very similar for Hill and MWC, while the optimal interacting networks are qualitatively different. 



We find that the general problem has many well de- 
fined local optima, corresponding to dilTerent topologies 
of the network, as shown in Fig |3] With just two target 
genes and a single input, the constraint of considering 
feed forward networks means that the locations of the ar- 
rows in the graph describing the network are fixed, but it 
is conventional to distinguish different topologies based 
on the signs associated to the arrows, that is whether a 
transcription factor acts as an activator (A) or repressor 
(i?), leaving us with multiple possibilities. For conve- 
nience we describe, for example, a network in which the 



input activates both target genes, while gene 1 represses 
gene 2, as A1A2 — Ri2- 

Our intuition is that the role of interactions is to re- 
duce redundancy, and this is borne out by the solutions 
associated to each topology. In particular, if the input 
activates both targets, then having one target activate 
the other can't help to reduce redundancy, and corre- 
spondingly we find that in the A1A2 — A12 topology the 
interaction between the target genes is actually driven to 
zero by the optimization, and this is true whether we use 
the Hill model of regulation or the MWC model (top and 
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bottom left panels of Fig|3|. The same intuition suggests 
that the R1R2 — A12 topology should also be driven to 
zero interaction, and this too is what we find (right panels 
of Fig [3]). As an aside, these comparisons also verify that, 
in the non-interacting case, the Hill and MWC models 
achieve very similar input/output relations once we find 
the parameters that optimize information transmission. 

The networks which are optimized with non~trival in- 
teractions are the ones where one target gene represses 
the other, A1A2 - R12 and RiAfi - R12. For the Hill 
model, both of these networks are optimized when one 
of the target genes is expressed only in a finite band or 
stripe of input concentrations, so that if we monitor only 
the expression level (72, we see a non-monotonic depen- 
dence on the input transcription factor. If we think of the 
input concentrations as varying along one axis of a two 
dimensional space — as they do in the genetic networks 
controlling embryonic development [T5] HTJ [TH] H^] — then 
this non-monotonic behavior means that one of the genes 
will be expressed in a "stripe" through the middle of the 
spatial domain. The other gene, which is itself a repres- 
sor, is expressed only at the extremes of the concentra- 
tion range, either at high or low input concentrations 
depending on whether the input is an activator or re- 
pressor, respectively. As in the case of non-interacting 
target genes, all of this structure depends on the available 
dynamic range of inputs. As we make Cmax/co smaller, 
the optimal strength of the interactions becomes smaller, 
and the networks that optimize information transmission 
are more nearly non-interacting; if Cmax/co is sufficiently 
small, the optimal solution again is for all targets to 
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FIG. 4: Optimal information transmission in different net- 
work topologies, with Hill model regulation. The information 
capacity of a network depends on Z, from Eq (381, which in 



turn is proportional to the number of independent copies of 
the output molecules, Ng. Here we plot Z — 2TTeZ/Ng as a 
function of the maximal concentration of input transcription 
factor, Cmax/Co. 
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FIG. 5: One target or two? We show the dependence of Z on 
the maximum concentration of input molecules, in the case of 
Hill regulation functions, comparing an interacting two gene 
network, a non-interacting two gene network, and a single 
output gene with twice the number of molecules. Note that 
to make the comparison meaningful, we have to be careful 
about the definition of co, since this concentration scale in- 
cludes a factor of Ng ; here we normalize to co for the two gene 
networks. 



be completely redundant, allowing the system to make 
multiple independent measurements of the same signal, 
rather than using the different targets to respond to dif- 
ferent portions of the input dynamic range. 

We know from the analysis of the non-interacting sys- 
tem that sub-optimal settings of all the parameters really 
do incur large penalties, substantially reducing informa- 
tion transmission [15j . This remains true for the inter- 
acting case, in that if we choose a topology but set the 
parameters arbitrarily, there is a significant loss of infor- 
mation. On the other hand, once we choose parameters 
optimally, the different topologies have very similar in- 
formation capacities, at least in the case of two target 
genes, as shown in Fig|4j To set a scale for these results 
on information transmission, we could ask whether there 
really is an advantage to having two target genes. After 
all, noise is reduced by making more molecules, and with 
two targets the maximal number of output molecules will 
be increased. A useful comparison, then, is with a system 
that has just one target gene, but can generate twice as 
many output molecules, and this is shown in Fig [5] We 
see that, as the structure of the optimal networks col- 
lapses at low values of c^h^/cq, so too does the advan- 
tage of having two separate target genes. At very low 
concentrations of the input transcription factors, it re- 
ally is better to have one target with twice the number 
of molecules, but as soon as there are enough molecules 
available to favor two target genes, the interactions be- 
tween these genes become important in building the net- 
works that optimize information flow. 
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A1A2AS — -R12-R13-R23 ^1^2^3 — -Rl2-Rl3^23 — ^12^13-^23 ^1^2^3 " ^12^13^23 
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FIG. 6: Locally opti ma l networks with the three target genes, with the input acting as an activator. Regulation is according to 
the Hill model, Eq (19 1, and we choose parameters to maximize information transmission with Cmax/co — 10. Top row shows 
the network topologies, and the middle row shows the average expression levels of the different target genes, as in Fig ID The 
bottom panel shows the distribution of input concentrations that is matches the optimal networks, from Eq (|37|). 



One of the most interesting results of the optimiza- 
tion problem is that the Hill and MWC models, which 
are almost indistinguishable in the non-interacting case, 
are driven to qualitatively different solutions in the in- 
teracting case. The basic idea of how interactions relieve 
redundancy can be understood in a limit where each gene 
is either fully expressed (5 = 1) or completely turned off 
{g = 0). For the non-interacting network, for example 
with the A1A2 topology, the genes "turn on" in sequence, 
and so the system accesses states 00, 01, and 11 in order 
as the input concentration is increased. In the interact- 
ing network A1A2 — R12, the optimal Hill model (cf Fig 
[3] top left panel) accesses state 00, 01, and then switches 
to 10, not quite achieving the state 11. In contrast, the 
MWC model, once optimized, accesses all the states in se- 
quence, 00 ^ 01 ^ 10 ^ 11, as the input concentration 
is increased, as we can see from the lower left panel of Fig 
[3] the situation is reversed in the case of R1R2 — Ri2- We 
can understand how this happens, because in the MWC 
model each binding event shifts the equilibrium between 
the active and inactive states; thus there is a competi- 
tion on gene 2 between the activating effect of the input 
and the repressing effect of gene 1, and this competition 
continues throughout the range of input concentration. 
In the MWC model, by adjusting parameters, it is pos- 
sible that the final increase in input concentration up to 



Cmax overwhelms the repressive effect of gene 1, and in- 
deed this is what happens at the parameter values that 
optimize information transmission. 

The difference between the Hill and MWC models indi- 
cates that the relatively microscopic physics of protein- 
DNA interactions that underlies transcriptional regula- 
tion can be "felt" at a more macroscopic, phenomenolog- 
ical level. In this spirit, we should also ask if it matters 
that we assume each protein can act only as an activator 
or a repressor, but not both. It is easy to imagine the 
alternative in which, by proper placement of the binding 
sites, the transcription factor can activate one gene by 
recruiting the polymerase, but can repress other genes 
by occlusion. We have analyzed these mixed networks, 
and find that the optimal information transmission in 
these cases is intermediate between the "pure" networks 
considered thus far. Thus, considering these additional 
topologies does not change the our global picture of the 
problem, except to open the possibility of yet more local 
optima. 



B. Three target genes 

The analysis of networks with three target genes con- 
firms and amplifies the lessons learned in the two gene 
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R1R2R3 — R12R13R23 R1R2R3 — -Rl2-Rl3^23 R1R2R3 — ^12^13-^23 
C C C 



R1R2R3 — ^12^13^23 




FIG. 7: Locally optimal networks with the three target genes, with the input acting as a repressor; all else as in Fig [6] 



case. Because the space of possible networks is much 
larger, we once again restrict individual proteins to act 
only as activators or repressors, and give explicit results 
only for the case of the Hill function models of regulation, 
Eq (19); results are shown in Figs |6] and [t] As with two 
genes, having the target genes activate one another can- 
not reduce the redundancy that is generated in the non- 
interacting networks, and hence the A1A2A2—A12A13A23 
topology networks are driven back to the non-interacting 
A1A2A3 when we solve the optimization problem, and 
similarly for R1R2R3 — ^12^13^23 Topologies that in- 
clude some repressive interactions generate non-trivial 
solutions, with progressively richer structure the more 
repression we allow. 

The signature of the repressive interactions in the case 
of two targets was the emergence of non-monotonic de- 
pendencies of the expression levels on the input transcrip- 
tion factor concentration. In the case of three targets, 
the parameters which maximize information transmission 
again lead to non-monotonicity. In the case where the 
input is an activator (Fig [6]), allowing for one repressive 
interaction {AiA2A3~Ai2Ai3R23) leads to one gene hav- 
ing a non-monotonic response. If we allow two repressive 
interactions {A1A2A3 — R12R13A23) then two genes have 
non-monotonic responses. Finally, with the maximum 
of three repressive interactions {A1A2A3 — R12R13R23) 
one of the two non-monotonic responses becomes yet 
more complex, with two "stripes" of expression in dif- 



ferent ranges of the input concentration. In the case of 
an input repressor, we see a similar pattern (Figj?]). 

A striking feature of the networks with two target 
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FIG. 8: Dependence of Z on the maximum available protein 
concentrations Cmax/co for two and three gene networks. We 
show results for the non-interacting cases, and for the globally 
optimal topologies, using both the Hill and MWC models of 
regulation. 
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FIG. 9: Optimal parameters for Hill regulatory functions in the globally optimal topologies of two {A1A2 — R12 grey lines) and 
three {A1A2A3 — R12R13R23 black lines) gene networks compared to the parameters of the optimal noninteracting networks 
{A1A2, A1A2A3). Parameters are presented as a function of the maximum concentration of available input proteins Cmax/co- 



genes was that local optima with different topologies 
were nearly degenerate. This is less true for the case 
of three targets — the spread in capacities associated 
with the different topologies is larger, and the network 
A1A2A3 — R12R13R23 is the clear global optimum. Cor- 
responding to the larger spread in capacities, the en- 
hancement of information transmission by interactions 
is larger, and in the difference in capacity between the 
Hill and MWC models also is larger in the three gene 
networks (Fig|8]). Including the possibility that individ- 
ual transcription factors could act as both activators and 
repressors does not change these conclusions. 

As a caveat, we note that we have assigned the same 
maximal concentration to all the transcription factors in 
the network. This is equivalent to assuming that all these 
molecules come at the same cost to the organism. One 
could imagine that there are significant differences be- 
tween the molecules, and hence that one could be ex- 
pressed at higher levels than the others for the same cost. 
If this symmetry is broken explicitly, then the ordering of 
the local optima associated to different network topolo- 
gies can be shifted. While it is satisfying to identify a 
global optimum, it probably thus still makes sense to 
think of the problem as having many local optima that 
might all be relevant in different biological contexts. An 
interesting feature of these different optima is that, de- 
spite the changing structure of the input / output relations 
for the target genes, the distribution of input concentra- 
tions that optimizes information transmission [cf Eq (37l] 



is almost the same for all the cases where the input is an 
activator (bottom panels of Fig [6]) and again for all the 
cases where the input is a repressor (bottom panels of 
Fig|7]). Although maximizing information flow requires 
matching of the input distribution to the characteristics 
of the regulatory network 6, 7 , it seems that once the 



parameters of the network itself have been optimized, the 
solution to the matching problem has a more universal 
structure, with almost the same distribution providing 
the best match to several different local optima. 

Finally, we consider the parameter values in these op- 
timal networks. For simplicity we focus on the globally 
optimal solutions with Hill model regulation functions for 
two {A1A2 — R12) and three {A1A2 — R12R13R23) tar- 
get genes, plotting in Fig [9] the evolution of the different 
constants Ki and Kij as a function of the maximal input 
concentration Cmax/co. As noted above, the optimal the 
interacting networks collapse onto the non-interacting 
case at small Cmax! for repressors, interactions become 
negligible if Kij becomes large, and this is what we see 
at small Cmax- At the opposite extreme, as Cmax/co be- 
comes large, the optimal Ki become nearly constant frac- 
tions of Cmax, while the interactions become stronger and 
stronger even up to Cmax/co ~ 10^. Put another way, the 
concentrations at which the target genes are turned on 
by the input become distributed at fixed fractions of the 
available dynamic range, while the points at which the 
targets are turned off by their mutual repression evolves 
as a function of the maximal input concentration. Fur- 
ther, in the interacting networks, the difference between 
two and three genes is not simply an extra target, but a 
readjustment of the interactions to make optimal use of 
this extra output channel. 



IV. DISCUSSION 

The expression levels of target genes provide informa- 
tion about the concentration of the transcription factors 
which provide input to the regulatory network. Our 
goal here has been to understand how the structure of 
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the network can be adjusted to optimize this informa- 
tion transmission. Earlier work [I5j addressed this prob- 
lem in the hmit where the target genes do not interact, 
so the input is simply "broadcast" to multiple targets; 
here we considered the role of interactions within a net- 
work of target genes, limiting ourselves, for simplicity, 
to feed forward structures. We have found that there 
are an exponentially large number of locally optimal net- 
works, each with parameters tuned to a well defined point 
that depends only on the maximal number of available 
molecules. These optimal networks have two compo- 
nents. As in the non-interacting case, the affinities of the 
input transcription factor for its multiple targets must 
be adjusted to balance the need to use the full dynamic 
range of inputs against the need to avoid noise at low in- 
put concentrations and output expression levels. Unique 
to the interacting case, mutual repression among the tar- 
get genes is tuned to reduce the redundancy among the 
network's outputs. 

The problem of redundancy reduction has a long his- 
tory in the context of neural coding, dating back to the 
first decade after Shannon's original work [51ji52j^S]. As 
first emphasized for the retina, nearby neurons often re- 
ceive correlated signals; since the dynamic range of neural 
outputs is limited, this redundancy compromises the abil- 
ity of the system to transmit information. The solution 
to this problem is lateral inhibition — neighboring neurons 
that receive correlated inputs inhibit one another, so that 
each neuron transmits something which approximates the 
difference between its input and that of its neighbors. 
Lateral inhibition, sometimes mediated through a more 
complex network, is a common motif in neural circuitry, 
and in the retina it is expressed by the "center-surround" 
organization of the individual neurons' receptive fields. A 
more careful analysis shows that there is a tradeoff be- 
tween redundancy reduction and noise reduction, so that 
the extent of inhibitory surrounds should be reduced, 
even to zero, as background light intensity is reduced 
and photon shot noise is increased, and this is in semi- 
quantitative agreement with experiment |54l 155] . The 
relation of these ideas in neural coding to our present 
discussion of genetic networks should be clear. Although 
there are many questions about how these ideas con- 
nect more quantitatively to experiment, it is attractive 
to think that similar physical principles, and even analo- 
gous implementations of these principles, might be rele- 
vant across such a wide range of biological organization. 

There are relatively few transcriptional regulatory net- 
works that have been characterized to the point of com- 
paring quantitatively with the sorts of models we have ex- 
plored here. In contrast, there is a growing literature on 
the qualitative, topological structure of these networks. 
In particular, a number of groups have focused on the 
local "motif" structure of large networks, searching for 
patterns of interaction that are over-represented relative 
to some randomized ensemble of networks |56l ISl| • One 
such motif is the feed forward loop, which essentially cap- 
tures the whole set of feed forward networks with two tar- 



get genes that we have considered here. In the language 
of Fig[l]D, much of the discussion about the importance 
of this motif has centered on the nature of the transfor- 
mation from the input concentration c to the expression 
level of gene 2. By choosing coherent (in our notation, 
AiA2~ A12) or incoherent {A1A2 — R12) loops, the system 
can achieve different temporal dynamics as well as non- 
monotonic input/output relations [58^,^, 60, 61, 62, 63]. 
Since we have considered system in steady state, it is this 
last point which is significant for our discussion. Our 
results comparing Hill and MWC models of regulation 
emphasize, however, that it is not just the topology of 
the network which is important for the structure of these 
input/output relations. 

One of the best opportunities for quantitative compari- 
son with experiment is in the genetic networks controlling 
the early events of embryonic development. In these sys- 
tems, the primary morphogen molecules are sometimes 
transcription factors, so that we can literally see the vari- 
ations in input concentration laid out in space [TTJ UHl ES] • 
Along the anterior-posterior axis of the Drosophila em- 
bryo, for example, the maternally provided morphogen 
Bicoid varies in concentration almost exponentially, so 
that the long axis of the embryo is essentially a logarith- 
mic concentration axis [SH [BS] . Bicoid is a transcription 
factor that provides input to the network of gap genes, 
hunchback, kriippel, giant and knirps. Although Bicoid 
is an activator, these genes exhibit non-monotonic ex- 
pression levels, forming stripes along the length of the 
embryo, and these non-monotonicities are mediated by 
mutual inhibitory interactions |57J [SHI IMl EO] . Qual- 
itatively, the structure of the optimal networks that we 
derive here (especially the globally optimal A1A2A3 — 
R12R13R23) resembles the gap gene network, and the ex- 
pression profiles along the log c axis are reminiscent of 
the spatial profiles of gap expression. Generalizations of 
the experiments in Ref [IT should make it possible to 
map the input / output relations more quantitatively, and 
perhaps even to measure the interactions directly by de- 
tecting the predicted correlations among the noise in the 
expression levels of different gap genes. To make truly 
quantitative comparisons, however, we need to solve our 
optimization problem allowing for networks with feed- 
back. 



APPENDIX 

This Appendix collects some technical details in the 
calculation of the inverse covariance matrix /C, following 
the framework established in Section |IID[ We start by 
considering the case of one input transcription factor at 
concentration c controlling the expression of two output 
genes, with normalized expression levels gi and g2- In 
this case the Langevin equations of motions describing 
the network are: 

^% = /i(c)-5i + 6(t), (53) 



dg2 



= /2(c,5l)-.92+C2(t)- 



(54) 



We linearize the equations of motion around the steady 
stat e so lutions {51,52} and Fourier transform to obtain 
Eq (46 1, where the matrix A is explicitly given by: 
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and hence 
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(55) The covariance matrix then can found, from Eq (52 1 
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The elements of A'ij are given explicitly by: 
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where we have chosen the maximum concentrations for 
all types of proteins to be equal to Cmax, and we have set 



Co = 1, as described in the text. To finish the calculation 
we need the integrals 
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Then 
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Lastly, we substitute the obtained covaiance matrix into 
the expression for the normalization constant Z, as ob- 



tained in Eq ( 38 1 : 
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(66) 



Since the optimal information is proportional to Z , we 
now can optimize Z over the parameters of the regula- 
tion function. It is now clear, having set the concentra- 
tion variables in terms of the natural scale cq = 1, that 



the only parameter in the problem is c^ax; all other pa- 
rameters will be determined by maximizing Z . 

In the case of a three gene network, the calculation is 
analogous to the two gene case. The matrix A is explicitly 
given by: 



A{u) = 



l-iujT 

— (/)21 1 — iujT 
-(f>31 -021 1 - il^T 



(67) 



and the elements of the noise spectrum matrix are 



16 



33 



1 



dgijc) 

dc 



c + gi 



df 2(0,91) 
d 13(0,91,92) 



c+ 021 1- 92 



dc 



,2 91 
hi- — 



{9.=9.(c)}, 



,2 92 
*32- 



'53 



(68) 
(69) 

(70) 



Following Eq (52 1 to find the inverse covariance integral, analogously to the two gene case, requires evaluating the 

(71) 



frequency integrals. The general form of these integrals is: 

duj 1 



'fcs 



2tt (1 + itJr)'^(l — ilot)'' 
1 f dC 1 

T 



1 /s + fc-2 



27r (1 + 0'=(1 - 



■yl — S — k 



TV fc- 1 

When the dust settles, the elements of the inverse of the covariance matrix, K,^^ from Eq (52 1, are: 
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^31 
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k — 1 
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^33^ 
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03l)'] + 


3 
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(72) 
(73) 

(74) 

(75) 

(76) 

(77) 

(78) 

(79) 

(80) 
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